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Abstract 

The equation modelling the evolution of a foam (a complex porous medium consisting 
of a set of gas bubbles surrounded by liquid films) is solved numerically. This 
model is described by the reaction-diffusion differential equation with a free boundary. 
Two numerical methods, namely the fixed-point and the averaging in time and 
forward differences in space (the Crank-Nicolson scheme), both in combination with 
Newton’s method, are proposed for solving the governing equations. The solution of 
Burgers’ equation is considered as a special case. We present the Crank-Nicolson 
scheme combined with Newton's method for the reaction-diffusion differential equation 
appearing in a foam breaking phenomenon. 

2000 Mathematics subject classification: primary 35K55; secondary 35R35, 65M06, 
65M12, 76B45, 76D99. 

Keywords and phrases: foam drainage, nonlinear partial differential equation, Crank- 
Nicolson method, breaking front, Plateau border, reaction-diffusion problem, free 
boundary problem. 


1. Introduction 

Foam drainage is a complex physico-chemical hydrodynamic process, as defined and 
explored in Kruglyakov et al. [16]. It is an important application area in industry since 
many technological properties of foam depend on the liquid volume fraction and the 
drainage rate. It is therefore important to understand how a foam forms and how it 
evolves in time. 

Drainage is defined as the transport of liquid through foam driven by gravity or 
pressure differences. The Plateau border (named after a blind Belgian scientist), which 
is a transition zone separating bulk surfaces, frames or other films from the him proper 
and always containing some bulk liquid, plays an important role in foam drainage 
studies. 
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Three factors, namely, the liquid distribution between the Plateau borders, the liquid 
films between two foam bubbles, and the distribution along the foam height, have to 
be known for the design of some technological processes; see, for example, Banhart 
and Weaire [2], Exerowa and Kruglyakov [11] and Weaire [19]. 

In Exerowa and Krugly akov [11], foam and foam films are considered, giving 
special attention to stability and drainage of foams. Experimental and theoretical 
aspects of the foam drainage equation are also considered in Koehler et al. [13], where 
an energy argument takes into account viscous dissipation in both the Plateau border 
and nodes to generalize the foam drainage equation. Three types of experiments, 
namely forced, free and pulsed drainage, are investigated. 

The numerical solution of the foam drainage equation has been studied by 
many researchers, including Colin and Fabrie [4, 5], Cox et al. [6], Durand and 
Langevin [10], Koehler et al. [14], Verbist et al. [18] and Weaire et al. [22]. In 2003, 
Weaire et al. [22] looked at problems concerning the fluid dynamics of foam under 
drainage, that is, the transport of liquid through foam. That the flow of liquid through a 
foam occurs mainly in its Plateau borders is shown using the surface evolver program 
(see Brakke [3]). An expanded review focusing on recent works on foam drainage 
including both the advanced theoretical and experimental studies of the foam drainage 
equation is given in Kruglyakov et al. [16]. Recently, Dahmani et al. [8] solved 
the foam drainage equation with time-fractional and space-fractional derivatives by 
applying the Adomian decomposition method. 

Here, however, we are interested in the case where the foam breaks due to the 
rupture of thin films during the drying process (see Colin and Fabrie [5], where 
the theoretical aspects of the convergence and well-posedness of the foam drainage 
equation are discussed). 


2. Formulation of the problem 

The foam drainage equation is a nonlinear partial differential equation that describes 
the evolution of the vertical density profile of foam under gravity. It was derived in 
1965 by Leonard and Lemlich [17] on the basis of various simplifying assumptions 
regarding the nature of this drainage process. The attribution of drainage entirely to 
Poiseuille flow in the network of Plateau borders is very important. These are the 
liquid-filled interstices of the bubbles which pack together to form the foam. 

Most of the essentials of the foam theory are to be found in Leonard and 
Lemlich [17]. In Gol’dfarb et al. [12], the following foam drainage equation is 
derived from the equation of continuity in the dimensional form, treating the liquid 
as incompressible for a single vertical Plateau border (channel-dominated) with cross- 
section A(r, §), which depends on the downward vertical coordinate £ and time r, 

9 9 

—A(r,£) + — [A(T,^(r,§)] = 0 (2.1) 

3r 

where the velocity q is averaged over the cross-section of the Plateau border. 
Gol’dfarb et al. [12] assumed that the flow in the channels formed by the Plateau 
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borders is of the Poiseuille type, that is, with zero velocity at the boundaries. 
In particular, the assumption of Poiseuille flow requires a sufficient surface viscosity to 
render the surfaces of the Plateau borders effectively rigid; see Kraynik [15] for details. 
Disregarding the film flow may also depend, to some extent, on this assumption; 
see Verbist et al. [18]. We also assume that drainage is sufficiently rapid that the 
diffusion of gas between cells is negligible. This is the process that leads to foam 
coarsening over long times. We neglect inertial effects and assume that there is a 
simple Poiseuille-type dissipation which is proportional to the mean liquid velocity in 
the Plateau border and inversely proportional to the cross section area. 

The following foam drainage equation is derived from assumptions which best 
apply to dry foam. The contribution of the films to drainage is not included in this 
model. On the basis of simple considerations of continuity and pressure balance, as 
well as the treatment of dissipation in a manner analogous to Darcy’s law in the theory 
of porous media, we follow the steps given in Weaire and Hutzler [21] to write the 
above Equation (2. 1 ) in dimensionless variables. 


3 u 
37 


+ 



= 0 


( 2 . 2 ) 


where 

vM 

2itpgy/nAi 

and where M = (>/3 — 7t/2) 1//2 is a geometrical constant corresponding to the 
triangular form of the Plateau border cross section and A,- is a typical value of the 
area of the Plateau borders. The physical parameters are v (constant surface tension), 
not including any effects of surface elasticity, p (liquid density), g (acceleration due 
to gravity), and 77 * (effective viscosity). (For more details see the appendix in Verbist 
etal. [18].) In Equation (2.2), u (t, x ) is the dimensionless cross section area of Plateau 
borders at time t and level x. that is, u(t, x) is a local measure of the foam’s liquid. 
Here, dimensionless coordinates, corresponding to vertical position x (measured 
downwards) and time t, are introduced by defining x = %/xq and t = t/ to where the 
units are xq = (Cv/pg) 1 ^ 2 and to = rj*/(Cvpg) 1 ' 2 , and u(t, x) = A(t, x)/xf y 

Note that Xq is (to within a constant) the capillary constant which has commonly 
been used in the theory of surfaces, and the height of meniscus is a useful indication 
of the magnitude of xq. The constant to, which sets the time-scale, is less familiar but 
is proportional to viscosity. 

The solution of (2.2) describes (for fixed t) a snapshot of the vertical density profile 
of the foam. The term in brackets in (2.2) is the dimensionless flow rate 


Q = u 2 



For a foam in equilibrium, we set <2=0, which gives an equilibrium profile 


Ueq(x) 


1 + R ~ X \ 

2 e ) 


-2 
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where u r is the liquid fraction at the bottom of the foam (that of a bubbly liquid at the 


rigidity loss transition where the bubbles of the foam separate) and R is the height of 
the foam. Notice that ur is never zero. 

A trivial solution of (2.2) is that of steady drainage with constant flow rate. Here 


it = ho, a constant, and so the flow rate is Hq and the flow velocity is hq. 


Equation (2.2) was called the foam drainage equation by Verbist et al. [18] and has 
to be completed with boundary conditions. In the case where the foam is stable during 
the draining (that is, does not break) we take Q = 0 at x = 0, and thus 



,2 


(2.3) 


u 


which means that no water enters the foam at the top level. The boundary condition 
at the bottom is more difficult to define. One reason for this is the assumption of 
dry foam while in reality the foam is wet at its bottom extremity where it touches the 
liquid. Here, we accept the assumptions that were suggested in Colin and Fabrie [5], 
where the problem posed on the half line was considered. Equation (2.2), associated 
with the above boundary condition (2.3), was considered in Koehler et al. [14]; we do 
not present any mathematical study of it, but draw attention to Koehler et al. [14] for 
numerical and experimental results. 

The situation where the foam breaks during the drying process is very important. 
In other words, at time t, the foam is confined in the domain \y(t). +oo), where y (t) 
denotes the position of the breaking front. In a physical sense, this means that the foam 
breaks when the area of the Plateau borders reaches a critical value > 0. Therefore, 
the boundary condition at x = y(t) is u(t, y(t)) = /?. 

So far, the case where dy(t)/dt = 0 has been examined. See Colin and Fabrie [5] 
for details, and the references therein. 

We consider the case where dy(t)/dt^0, which is extremely important and 
completely different from the cases considered up to now. It is known that conservation 
of the quantity of water at the breaking front leads us to 


dy(t) s f 3 u(t, x ) 


(2.4) 


yffi - dx _ x=y(t) 


Therefore, the problem we have to deal with reads 



(2.5) 


with the boundary conditions 
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and the initial condition 


m(0 , x) = uo(x), x € (y(f), +oo) 


where 0 < ft < 1 and u o(x) > ft, x e (y(0), +oo). In order to study (2.5) we make a 
transformation v(t, x) = u(t, x + y (f)). Then the equation satisfied by v is 

dv dv dv 2 3 ( / -dv\ 

-y(0—t-e— Vv— =o, 

3 1 r dx dx dx\ 3 xj 

with the boundary conditions 
v(t, 0) = p. 


m = f ~TpL 


3n(f, x) 
dx 


y(0) = o, 


Jx=0 


and the initial conditions 


y(0) =0, 

v(0, x) = vo(x), i e (0, +oo). 

We consider the bounded interval (0, R) with homogeneous Neumann boundary 
condition at x = R. Hence, the problem becomes 


dv dv dv 2 3 ( ,-dv\ 

- y(t) -1- e —( s/v — 1=0, t > 0, 

dt r dx dx dx\ dx) 


( 2 . 6 ) 


with the boundary conditions 
v(t, 0) = p. 


3n(r, x) 
dx 


= 0, t> 0, 


y( ° P dx 


J x=R 
dv(t, x) 


y (0) = 0, t> 0 


J^=o 


and the initial conditions 


y (0) = o, 

u(0, x) = i>o(x), x € (0, R) 


where R is the height of the foam and an arbitrary but fixed positive number. Colin 
and Fabrie [5] considered Equation (2.6) and assumed that 

0<p<\, 

vo e W 1,oo (0, R), vo(0) = P, 


vo(x) > P, x e (0, R) 

where W l ’°° is the usual Sobolev space; see Adams [1]. Colin and Fabrie [5] 
constructed the global solutions of (2.6). 
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We consider the equation governing the evolution of a foam Equation (2.6) with the 
above boundary and initial conditions. Cox et al. [7] considered the case of y(t ) = 0 
using a constant-step-size finite-difference representation, with forward differencing 
in time (upwind) and space. 


3. Numerical procedure 


We assume that v{t . x) > 0, given by (2.6), since v(t, x ) represents the 
dimensionless cross section area of Plateau borders at time t and level x. 

Therefore, we use the transformation z(f, x) = \ v(t, x) 3/2 ]. Hence, Equation (2.6) 
becomes 


2 1 dz 2 1 dzf 2s y dz(t,x) 

3x J x=0 

with the boundary conditions 

3 z(t, x) 

. dx J x=R 

2 1 3 z(t, x) 
3z 1//3 dx 




z(t, 0) = 


= 0, t > 0, 




Jx=0 


and initial conditions 


(3.2) 


t > 0 


y(0) = 0, t> 0, 
z(0, x) = zo(x), x G (0, R). 
Equation (3.1) is multiplied by ^z 1 / 3 , giving 

3 zit, x)~| \ „ on dz 


dz dz / 2 s 

8t ~ dx V ~~ 30 . 


dx 


d 2 Z 


+ 2z 2 ' 3 ^-sz 1 ' 3 ^ = 0. 
j x= q/ dx dx 


(3.3) 


3.1. Fixed-point iteration method The partial differential equation (3.3) is 
discretized in space by standard finite differences on an equally spaced grid and using 
a forward-time difference scheme with constant step size. The implicit scheme then 
becomes 


7 n +1 _ 7 n 
'■m _-■/// 

At 


7 n +1 _ n +1 
4 m+1 


Ax 


2g / z" +l — Zq + 1 
P 3/1 \ Ax 


2(z" +1 ) 2/3 ( 


n +1 


Z m+\ 

Ax 


n +1 


/ 7 » + l _ 2 ? n + l 4- 7 n + 1 
„/„«+l\l/3 / ^m+l ^'■hi- 1 

Um ; V (a^) 2 


= 0 


where z” = z(t«, •*»/) = z(nk , m/z), £ = At, h = Ax and w = 0, . 
refer to this method as Method I. 


(3.4) 
M - 1. We 
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Equation (3.4) becomes 


, 7 «+l _ n+l 
7 n+l _ 7 n , At l 4 m+l 
4 ra 


Ax 


/ 7 ' ,+ 1 — 7 ,!+1 

ft _ _/ ■“! 0 


3/6 \ Ax 


/z" 1 " 

-2Af(z ) ; +1 ) 2/3 f : ^ ±i 


7 n+l _ 7 «+l 


+ £ Af(z” +1 ) 1/3 (^ ±1 


Ax 

t"+ 1 _O+z+l _i_ 7 ' ! +l 

^m+1 ~ 1 


(Ax) 2 

Fixed point iteration is applied to (3.5) in time using 


n+l,k+\ _ n,k 


= 4’* + Ar 


n+l,k _ n+l,k 

Z m +1 


Ax / V' 3/6 
w+l,£ w+l,fc 




Ax 


-2AKz” +U ) 2/3 ( W ^ ' ) 

✓ w+l,£ ^ w+l,fc n+l,k\ 

+ e A >( C - t )'/ 3 ( Z - +1 ~ ( y j); ) 


which is written in more compact form as 


(3.5) 


n+l,fc+l _ n,k , cy n+l,* n+l,ft; n+l,*-. 
•c m 4 m T * —1 


This fixed point iteration method is first-order accurate. Thus, this implicit method 
(Method I) is only first-order accurate in time and second-order accurate in space. 

3.2. The Crank-Nicolson scheme The Crank-Nicolson scheme is applied to (3.3), 
yielding 


„n +1 _ 7 n 
'■in _ ''iii 

At 


7 »+l _ 7 n+ 1 
4 m+1 4 m 


Ax 


7 — 7 

•“Hl+1 


Ax 


7 « - 7 « \ 1 

41 4, 


- 2g l A” +1 — £q +1 <■, 

3/6 2 \ Ax Ax 


+ 


-71+1 


+ Z 

I +> ,1 


2/3 /_«+l _ n+l 
z m+ 1 7 


+ 


^m+1 


Zri 


£ / + 7^ 

^ I +• m 1 


2 y V Ax Ax 

n+i _ 27 ,!+1 + 7 n+l 7 n -2v n + 7 " 

'•m + 1 '■m ' "m -1 . '•ni+1 T<, m-1 


(Ax) 2 


+ 


(Ax) 2 


2 

= 0 . 


1/3 


We define the M mesh points x/= (/— l)£/(2 1 / 3 ), 7 = 1, ... , M where 
E = 0.11 and the interval of discretization is chosen arbitrarily. The +J+ 1 values are 
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TABLE 1. The range of the parameters for each simulation. 


Figure 

e 

P 

Ax 

At 

Method 

1 

0.01 

0.001 

0.0873 

0.0873 

Crank-Nicolson 

2 

0.01 

0.001 

0.0873 

0.1 

Implicit (Method I) 

5 

0.01 

0.99 

0.0873 

0.0873 

Crank-Nicolson 

6 

0.01 

0.99 

0.0873 

0.0873 

Implicit (Method I) 

7 

10~ 7 

0.001 

0.0873 

0.0873 

Crank-Nicolson 

8 

10” 7 

0.99 

0.0873 

0.0873 

Crank-Nicolson 


unknowns for m= 0, . . . , M — 1. For in = 0, we use the forward difference formula; 
for m = M — 1, we impose the boundary condition [ dz.it. x)/dx \ x _ R =0 using the 
backward second-order difference formula z n M = 4/3 z n M _ l — 1/3 z n M _~,- 

In both Method I and the Crank-Nicolson scheme for each time step we must solve 
M nonlinear algebraic equations to find M unknowns for fixed values of e and /l. We 
solve that system using Newton’s method. 

4. Discussion of results 

We use the numerical procedures described in Section 3 to compute the solutions 
for various values of /J and s. In most of the computations E = 0.11 and M = 400. 
We ran the calculations with different values of E and M and observed that the results 
are independent of E and M to graphical accuracy. 

In numerical calculations, we take all quantities in dimensionless form, so that we 
obtain the corresponding physical quantities by multiplying lengths by xq and time 
by t(). Table 1 shows the range of the parameters, namely e and /l, space and time 
increments Ax and At respectively, and the numerical methods. 

We take the uniform solution as an initial guess since it satisfies not only the given 
differential equation but also the boundary and initial conditions. 

The exact solution of the foam drainage equation has not yet been obtained, as far as 
we are aware. No tables were given by Cox et al. [7] or Weaire et al. [20, Page L219, 
Figure 6], but a comparison with their graphs shows good agreement. 

In our calculations, the convergence tolerance was 10 8 , that is, we terminated the 
iterative procedure when successive approximations agreed up to the eighth decimal 
place. 

In Figure 1, at the left-hand side (top of the foam) the liquid fraction tends to a 
constant, corresponding to a wet foam. At the bottom, the liquid fraction is almost 
zero, corresponding to a dry foam. 

We observe that the error in the Crank-Nicolson scheme has a quadratic 
convergence rate, and in Method I it has a linear convergence rate (see Figure 2). 
In Figure 1, it took 14000 iterations to reach the equilibrium solution. However, 
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FIGURE 1, u(t, *) versus* plotted using the Crank-Nicolson method: e = 0.01, ft = 0.001, Ax = 0.0873 
and At = 0.0873, for different times t = 5, 25, 55, 91.6, 218.2, 353.5 and 1222.3. The leftmost curve 
shows the solution at t = 5. As time progresses, the curves move from left to right. 



FIGURE 2. u{t, x) versus * plotted using Method I (implicit method): e = 0.01, ft = 0.001, Ax = 0.0873 
and At = 0.1, for different times t = 5, 25, 48, 130, 250, 436 and 1450. The leftmost curve shows the 
solution at t = 5. As time progresses, the curves move from left to right. 


in Figure 2 (using Method I), At = 0.0873 was too small to reach the equilibrium 
solution so At was taken as 0.1. The problem with Method I is that it is only first- 
order accurate in the time step and second-order accurate in space. On the other hand, 
the Crank-Nicolson method is second-order accurate in time and space. Notice that 
both methods are unconditionally stable. 
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FIGURE 3. u(t , x) versus x plotted using the Crank-Nicolson method: e = 0.01, fl = 0.99, At = 0.0873 
and Ax = 0.0873, for different times t = 0.0873 (at the end of the first step) and then t = 5.1, 7.7, 12.3, 
16.5, 21.9, 25.4, 34.1 and 42.8. The leftmost curve shows the solution at t = 0.0873. As time progresses, 
the curves move from left to right. 



X 


FIGURE 4. u(t, x) versus x plotted using Method I (implicit method): e = 0.01, ft = 0.99, Ax = 0.0873 
and At = 0.0873 for different times t = 0.0873 (at the end of the first step), 1.3, 2.7, 5, 7.6, 12.2, 16.4 
and 21.8. The leftmost curve shows the solution at t = 0.0873. As time progresses, the curves move from 
left to right. 


The foam drainage equations have a travelling wave solution. In Figures 3, 4 and 7 
the liquid fraction looks like a solitary wave, which is a wave of constant profile, 
moving through the foam with constant velocity. However, eventually all cases reach 
an equilibrium solution. 
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FIGURE 5. \y(t)\ versus t shown using the Crank-Nicolson method: e = 0.01, fi = 0.001, Ax = 0.0873 
and At = 0.0873. 



FIGURE 6. Zoomed version of \y(t)\ versus t shown using the Crank-Nicolson method: e = 0.01, 
P = 0.001, Ax = 0.0873 and At = 0.0873. 


4.1. The position of the breaking front: y ( t ) In Section 2 we described y (t ) as the 
position of the breaking front. The foam breaks when the area of the Plateau borders 
reaches a critical value fi > 0. The boundary condition at x = y(t) is u(t, y(t )) = fi. 
The function y(t) is computed from y (t), that is, 


y(t) = 


L 


r dy(s) 


ds. 


ds 
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FIGURE 7. u(t, x) versus x plotted using the Crank-Nicolson method: e = 10“ 7 , fl = 0.99, Ax = 0.0873 
and At = 0.0873 for different times t = 0.0873 (at the end of the first step), t = 4.4, 7, 12.3, 16.6, 21, 
25.4, 29.7 and 42.8. The leftmost curve shows the solution at t = 0.0873. As time progresses, the curves 
move from left to right. 



X 


Figure 8. u(t, x) versus x plotted using the Crank-Nicolson method: s = 10 1 , = 0.001, Ax = 

0.0873 and At = 0.0873 for different times t = 5, 25, 55, 87.3, 171.1, 433 and 1527.8. The leftmost 
curve shows the solution at t — 5. As time progresses, the curves move from left to right. 


As explained in Colin and Fabrie [5], the conservation of the quantity of water at the 
breaking front gives us the description of y(t), given in Equation (2.4). We calculate 
the above integral by applying the trapezoidal rule. 

The values of y (f) are positive and the curves of y(t ) behave like a linear function 
except for e = 0.01 and 0 < /3 < 0.2. Our calculation with (i in the range 0 < /3 < 0.2 
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shows that y (f) becomes negative. One reason might be that y (t) in Equation (3.2) is 
negative for small values of /?. The profiles corresponding to this case are shown in 
Figure 5 (the absolute value of y (f)) and enlarged in Figure 6. 

4.2. Special case: e = 10 -7 , Burgers’ equation Notice that the general foam 
drainage Equation (2.5) reduces to an equation similar to the well-known inviscid 
Burgers’ equation if s goes to zero; see Debnath [9], The two equations are not equal, 
since the inviscid Burgers’ equation has the form 

3 u 1 3 u 2 

37 + 2~37 = ’ 

including the factor 1/2. In Equation (2.6), s is taken to be 10 7 to compare with 
Burgers’ equation (see Figures 7 and 8). In both cases, y(t) behaves like a linear 
function. 


5. Conclusion 

We have solved numerically the equations of a one-dimensional model of the 
evolution of a foam using both an implicit method and the Crank-Nicolson method in 
combination with Newton’s method. The corresponding problem in two dimensions is 
left for future work. 
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